LMS, Signed-LMSの挙動観察

LMS, SIgned-LMSの挙動観察

入力信号を固定し、各フィルタ係数において、残差、残差勾配はどうなっているのか？ まずはいつもどおりの LMS フィルターからはじめ、Signed-LMS の解析を行ってみる。

In [1]:
import sys
import math
import numpy as np
import matplotlib.pyplot as plt

# 固定した係数による LMS フィルター処理
def lms_fixcoef(data, coef):
    num_samples = data.size
    filter_order = coef.size
    pred = np.zeros(num_samples)
    for smpl in range(filter_order, num_samples, 1):
        pred[smpl] = np.dot(coef, data[smpl - filter_order:smpl])
    return pred
In [2]:
# xtics, ytics で係数を生成して LMS フィルター処理を実行
# 各点における残差と残差勾配を出力
def lms_experience(data, xtics, ytics):
    zerror = np.zeros((xtics.size, ytics.size))
    zerrorgrad = np.zeros((xtics.size, ytics.size, 2))
    for x in range(xtics.size):
        for y in range(ytics.size):
            coef = np.array([xtics[x], ytics[y]])
            pred = lms_fixcoef(data, coef)
            error = data - pred
            zerror[x, y] = np.sqrt(np.mean(error ** 2))
            zerrorgrad[x, y, 0] = -2 * np.mean(error[1:] * data[0:-1])
            zerrorgrad[x, y, 1] = -2 * np.mean(error[2:] * data[0:-2])
    return zerror, zerrorgrad

def ploterror(label, xtics, ytics, zmesh):
    xmesh, ymesh = np.meshgrid(xtics, ytics)
    # plt.pcolormesh(xmesh, ymesh, zmesh, label=label, cmap='Blues')
    plt.contour(xmesh, ymesh, zmesh, levels=np.linspace(0.0, 1.5, 10).tolist())
    plt.colorbar()
    plt.xlabel('x0')
    plt.ylabel('x1')
    
def plotgradient(label, xtics, ytics, zmeshvec):
    xmesh, ymesh = np.meshgrid(xtics, ytics)
    plt.quiver(xmesh, ymesh, zmeshvec[:,:,0], zmeshvec[:,:,1], label=label)
    plt.xlabel('x0')
    plt.ylabel('x1')

# 実験用入力データの作成
NUM_SAMPLES = 1024
lineardata = np.ones(NUM_SAMPLES)
sindata = np.ones(NUM_SAMPLES)
for i in range(NUM_SAMPLES):
    sindata[i] = math.sin(2 * math.pi *  440 * i / 48000.0)

np.random.seed(0)
gaussnoise = np.random.normal(0, 0.1, NUM_SAMPLES)
np.random.seed(0)
laplacenoise = np.random.laplace(0, 0.1, NUM_SAMPLES)

# x,y の解析範囲
xtics = np.linspace(-2, 3, 20)
ytics = np.linspace(-2, 3, 20)

plt.plot(lineardata)
plt.show()
plt.plot(sindata)
plt.show()
plt.plot(gaussnoise)
plt.show()
plt.plot(laplacenoise)
plt.show()
plt.plot(sindata + gaussnoise)
plt.show()
plt.plot(sindata + laplacenoise)
plt.show()
In [3]:
zerror, zerrorgrad = lms_experience(lineardata, xtics, ytics)
        
ploterror('Error', xtics, ytics, zerror)
plotgradient('Error Gradient', xtics, ytics, zerrorgrad)
plt.title('Error and Error Gradient')
plt.show()
In [4]:
zerror, zerrorgrad = lms_experience(sindata, xtics, ytics)

ploterror('Error', xtics, ytics, zerror)
plotgradient('Error Gradient', xtics, ytics, zerrorgrad)
plt.title('Error and Error Gradient')
plt.show()
In [5]:
data = lineardata + gaussnoise
zerror, zerrorgrad = lms_experience(data, xtics, ytics)
        
ploterror('Error', xtics, ytics, zerror)
plotgradient('Error Gradient', xtics, ytics, zerrorgrad)
plt.title('Error and Error Gradient')
plt.show()
In [6]:
data = sindata + gaussnoise
zerror, zerrorgrad = lms_experience(data, xtics, ytics)

ploterror('Error', xtics, ytics, zerror)
plotgradient('Error Gradient', xtics, ytics, zerrorgrad)
plt.title('Error and Error Gradient')
plt.show()
In [7]:
data = lineardata + laplacenoise
zerror, zerrorgrad = lms_experience(data, xtics, ytics)
        
ploterror('Error', xtics, ytics, zerror)
plotgradient('Error Gradient', xtics, ytics, zerrorgrad)
plt.title('Error and Error Gradient')
plt.show()
In [8]:
data = sindata + laplacenoise
zerror, zerrorgrad = lms_experience(data, xtics, ytics)

ploterror('Error', xtics, ytics, zerror)
plotgradient('Error Gradient', xtics, ytics, zerrorgrad)
plt.title('Error and Error Gradient')
plt.show()

ここまでの考察

  • どんな波形を突っ込んでも同じ残差分布になる（直線 x0+x1=1 上に解があるように見える）...
    • 音声特有の相関の強みが出ている。ほぼ、x0+x1=1 が理想的になる。x1（2 個前データの係数）=1 とするよりも x0（直前データの係数）=1 の方が残差が小さい。これは、直前の信号により相関しているから。x0 = x1 = 0.5 はその平均になり、これも残差は小さい。
    • 勾配は正しく計算できてそう。
  • 雑音を付加すると残差分布が楕円状に変わる
    • 雑音が加わると、個々のサンプルの値よりも、サンプルの平均を取ったときに残差が小さくなる。
      • つまり、x0 = x1 = 0.5 が理想に近い。
    • 正規（ガウス）雑音、ラプラス分布雑音でほぼ同じ傾向。しかしラプラスの方がやや裾が広い？
    • ふつーの LMS はガウス雑音前提でやってるので注意。
    • 勾配の向きが怪しい。グラフ左上、右下が残差勾配を下っている（= 悪化している）様に見える。
    • → 実装不具合のようだ。予測時のフィルタ係数を内積と一致させる（coef[0] data[smpl - 2] + coef[1] data[smpl - 1] とさせ）ないと座標軸と一致しない（反転してしまう。）

→ 次はSigned-LMSで解析してみる。

In [9]:
# xtics, ytics で係数を生成して Signed-LMS フィルター処理を実行
# 各点における残差と残差勾配を出力
def signedlms_experience(data, xtics, ytics):
    zerror = np.zeros((xtics.size, ytics.size))
    zerrorgrad = np.zeros((xtics.size, ytics.size, 2))
    for x in range(xtics.size):
        for y in range(ytics.size):
            coef = np.array([xtics[x], ytics[y]])
            pred = lms_fixcoef(data, coef)
            error = data - pred
            zerror[x, y] = np.mean(np.abs(error))
            zerrorgrad[x, y, 0] = - np.mean(np.sign(error[1:]) * data[0:-1])
            zerrorgrad[x, y, 1] = - np.mean(np.sign(error[2:]) * data[0:-2])
    return zerror, zerrorgrad

data = lineardata
zerror, zerrorgrad = signedlms_experience(data, xtics, ytics)
        
ploterror('Error', xtics, ytics, zerror)
plotgradient('Error Gradient', xtics, ytics, zerrorgrad)
plt.title('Error and Error Gradient')
plt.show()
In [10]:
data = sindata
zerror, zerrorgrad = signedlms_experience(data, xtics, ytics)

ploterror('Error', xtics, ytics, zerror)
plotgradient('Error Gradient', xtics, ytics, zerrorgrad)
plt.title('Error and Error Gradient')
plt.show()
In [11]:
data = lineardata + gaussnoise
zerror, zerrorgrad = signedlms_experience(data, xtics, ytics)
        
ploterror('Error', xtics, ytics, zerror)
plotgradient('Error Gradient', xtics, ytics, zerrorgrad)
plt.title('Error and Error Gradient')
plt.show()
In [12]:
data = sindata + gaussnoise
zerror, zerrorgrad = signedlms_experience(data, xtics, ytics)

ploterror('Error', xtics, ytics, zerror)
plotgradient('Error Gradient', xtics, ytics, zerrorgrad)
plt.title('Error and Error Gradient')
plt.show()
In [13]:
data = lineardata + laplacenoise
zerror, zerrorgrad = signedlms_experience(data, xtics, ytics)
        
ploterror('Error', xtics, ytics, zerror)
plotgradient('Error Gradient', xtics, ytics, zerrorgrad)
plt.title('Error and Error Gradient')
plt.show()
In [14]:
data = sindata + laplacenoise
zerror, zerrorgrad = signedlms_experience(data, xtics, ytics)
        
ploterror('Error', xtics, ytics, zerror)
plotgradient('Error Gradient', xtics, ytics, zerrorgrad)
plt.title('Error and Error Gradient')
plt.show()

Signed-LMSの観察

  • 雑音がない時、どれだけ解から離れていても勾配が同一
    • 稜線が直線になった谷になっていると想像。
    • 普通の LMS では解から離れると 2 乗則？で勾配が大きくなる
    • 逆に、どんなに近くても勾配が同一（LMS では解の近傍で勾配が急激に消える）
    • 残差の裾野が LMS に比べて広い。
  • 雑音がある時でも、解の近傍で勾配が消失しにくい。
    • 解から離れると雑音がない時の勾配に近づく。
    • それでも尾根方向は勾配が消えている。
  • やっぱ山の稜線に沿った勾配（左上から右下）の向きが気になる。。。これ学習しようとすると思いっきり悪い方向に進むよな。
    • 勾配の向きが怪しいのは、勾配法では必ずしも解の稜線に沿った向きに動かないから？→ 実装バグだった。
    • ここ あたりが参考になるかも？勾配は必ずしも極値点を向いてない。

→引き続き、LMSの勾配に観測分散行列（ヘッセ行列）の逆行列を掛けたらどうなるか見ていく

In [15]:
# xtics, ytics で係数を生成して LMS フィルター処理を実行
# 各点における残差と残差勾配を出力
def lms_experience_hessian(data, xtics, ytics):
    zerror = np.zeros((xtics.size, ytics.size))
    zerrorgrad = np.zeros((xtics.size, ytics.size, 2))
    zerrorgradhessian = np.zeros((xtics.size, ytics.size, 2))
    H = np.zeros((2, 2))
    for smpl in range(2, data.size, 1):
        xvec = data[smpl - 2:smpl].reshape((2,1))
        H += xvec @ xvec.T
    H /= data.size
    Hinv = np.linalg.inv(H)
    for x in range(xtics.size):
        for y in range(ytics.size):
            coef = np.array([xtics[x], ytics[y]])
            pred = lms_fixcoef(data, coef)
            error = data - pred
            zerror[x, y] = np.sqrt(np.mean(error ** 2))
            grad = -2 * np.array([np.mean(error[1:] * data[0:-1]), np.mean(error[2:] * data[0:-2])])
            zerrorgrad[x, y, :] = grad
            zerrorgradhessian[x, y, :] = (Hinv @ grad.reshape((2,1))).reshape(2)
    return zerror, zerrorgrad, zerrorgradhessian, Hinv

xtics = np.linspace(-2, 3, 20)
ytics = np.linspace(-2, 3, 20)

# data = lineardata # 特異になる！
data = sindata # あやしい。なぜか (2,-1) 付近に集まっとる
zerror, zerrorgrad, zerrorgradhessian, Hinv = lms_experience_hessian(data, xtics, ytics)

ploterror('Error', xtics, ytics, zerror)
plotgradient('Hessian Gradient', xtics, ytics, zerrorgradhessian)
plt.title('Error and Hessian Error Gradient')
plt.show()
In [16]:
data = lineardata + gaussnoise
zerror, zerrorgrad, zerrorgradhessian, Hinv = lms_experience_hessian(data, xtics, ytics)

ploterror('Error', xtics, ytics, zerror)
plotgradient('Hessian Gradient', xtics, ytics, zerrorgradhessian)
plt.title('Error and Hessian Error Gradient')
plt.show()
In [17]:
data = sindata + gaussnoise
zerror, zerrorgrad, zerrorgradhessian, Hinv = lms_experience_hessian(data, xtics, ytics)

ploterror('Error', xtics, ytics, zerror)
plotgradient('Hessian Gradient', xtics, ytics, zerrorgradhessian)
plt.title('Error and Hessian Error Gradient')
plt.show()
In [18]:
data = lineardata + laplacenoise
zerror, zerrorgrad, zerrorgradhessian, Hinv = lms_experience_hessian(data, xtics, ytics)

ploterror('Error', xtics, ytics, zerror)
plotgradient('Hessian Gradient', xtics, ytics, zerrorgradhessian)
plt.title('Error and Hessian Error Gradient')
plt.show()
In [19]:
data = sindata + laplacenoise
zerror, zerrorgrad, zerrorgradhessian, Hinv = lms_experience_hessian(data, xtics, ytics)

ploterror('Error', xtics, ytics, zerror)
plotgradient('Hessian Gradient', xtics, ytics, zerrorgradhessian)
plt.title('Error and Hessian Error Gradient')
plt.show()

ヘッセ行列を使った時の観察

  • ノイズがない場合、逆行列が計算できない場合がある。
    • 分散行列が特異になる。見たところ分散行列の全ての要素が同一。
    • 理論的には確かに半正定値までしか保証されていないから、正則ではないのはありえる。
    • 雑音を入れた場合は今の所逆行列が計算できてる。でも、理論上は半正定値なので常に計算できるとは限らず。
    • サイン波ではノイズなしで逆行列を計算できたけど勾配が (2, -1) 付近を向いてる。
  • ノイズがある場合、全ての勾配が極値点 (0.5, 0.5) 付近を向いている。
  • 勾配の尾根方向における勾配消失が見られず、純粋に極値からの距離で勾配が決まっている印象がある。
    • これは最適化で有利に働きそう。

→引き続き、Signed-LMSに対して実験してみる

  • 定義式通りインパルス関数を入れたらどうなってしまうんだろうか …（緊張）
In [20]:
# 平均 0 分散 var のラプラス分布の P(|x| <= delta) を計算
def laplace_distribution_underdelta(var, delta):
    return 1 - math.exp(-delta/var)

# ラプラス分布の 0 とみなす閾値
DELTA = 0.01

# xtics, ytics で係数を生成して LMS フィルター処理を実行
# 各点における残差と残差勾配を出力
def signedlms_experience_hessian(data, xtics, ytics):
    zerror = np.zeros((xtics.size, ytics.size))
    zerrorgrad = np.zeros((xtics.size, ytics.size, 2))
    zerrorgradhessian = np.zeros((xtics.size, ytics.size, 2))
    for x in range(xtics.size):
        for y in range(ytics.size):
            coef = np.array([xtics[x], ytics[y]])
            pred = lms_fixcoef(data, coef)
            error = data - pred
            H = np.zeros((2, 2))
            count = 0
            for smpl in range(2, data.size, 1):
                xvec = data[smpl - 2:smpl].reshape((2, 1))
                H += xvec @ xvec.T
            var = np.mean(np.abs(error))
            H = H * laplace_distribution_underdelta(var, DELTA) / data.size
            Hinv = np.linalg.inv(H)
            zerror[x, y] = np.sqrt(np.mean(error ** 2))
            grad =  -np.array([np.mean(np.sign(error[1:]) * data[0:-1]), np.mean(np.sign(error[2:]) * data[0:-2])])
            zerrorgrad[x, y, :] = grad
            zerrorgradhessian[x, y, :] = (Hinv @ grad.reshape((2,1))).reshape(2)
    return zerror, zerrorgrad, zerrorgradhessian

# data = lineardata # 特異になる！
data = sindata
zerror, zerrorgrad, zerrorgradhessian = signedlms_experience_hessian(data, xtics, ytics)

ploterror('Error', xtics, ytics, zerror)
plotgradient('Hessian Gradient', xtics, ytics, zerrorgradhessian)
plt.title('Error and Hessian Error Gradient')
plt.show()
In [21]:
data = lineardata + gaussnoise
zerror, zerrorgrad, zerrorgradhessian = signedlms_experience_hessian(data, xtics, ytics)

ploterror('Error', xtics, ytics, zerror)
plotgradient('Hessian Gradient', xtics, ytics, zerrorgradhessian)
plt.title('Error and Hessian Error Gradient')
plt.show()
In [22]:
data = sindata + gaussnoise
zerror, zerrorgrad, zerrorgradhessian = signedlms_experience_hessian(data, xtics, ytics)

ploterror('Error', xtics, ytics, zerror)
plotgradient('Hessian Gradient', xtics, ytics, zerrorgradhessian)
plt.title('Error and Hessian Error Gradient')
plt.show()
In [23]:
data = lineardata + laplacenoise
zerror, zerrorgrad, zerrorgradhessian = signedlms_experience_hessian(data, xtics, ytics)

ploterror('Error', xtics, ytics, zerror)
plotgradient('Hessian Gradient', xtics, ytics, zerrorgradhessian)
plt.title('Error and Hessian Error Gradient')
plt.show()
In [24]:
data = sindata + laplacenoise
zerror, zerrorgrad, zerrorgradhessian = signedlms_experience_hessian(data, xtics, ytics)

ploterror('Error', xtics, ytics, zerror)
plotgradient('Hessian Gradient', xtics, ytics, zerrorgradhessian)
plt.title('Error and Hessian Error Gradient')
plt.show()

SignedLMS・ヘッセ行列を使った時の観察

怪しいけど、残差 0 のときに和を取るのをやめて、残差 0 を取る確率で重み付けして分散行列を計算した。

  • 尾根の平坦なところ（画面左上、右下）で大きい勾配が出るようになった。
  • 反面、左下と右上で勾配が小さくなったように見える（大きさは変わってないけど、相対的に小さくなった？）
  • 分散パラメータ 0.3 は小さすぎる気がする。（残差が 0 になる確率が 0.93 と大きすぎる。）離散のスケールにそぐわないような気がする。
    • はじめから離散スケールで実験を行うべきかも。信号の振幅は [-32768, 32767] とする。
    • いや、[-1,1] のスケールで、残差 0 確率が 10% になるように分散を設定するのが正しいか。
    • → 連続型確率分布で、残差の絶対値が閾値 delta 以下になる確率を計算するようにした

→実際に学習をシュミレーションしてみる

  • 初期値を色々変えたときにどこに行くか。
  • ヘッセ行列を使った時に最適値に達する回数は？
In [25]:
# LMS フィルター処理
def lms(data, initcoef, stepsize):
    num_samples = data.size
    filter_order = initcoef.size
    pred = np.zeros(num_samples)
    error = np.zeros(num_samples)
    coef = initcoef
    coefhistory = np.zeros((filter_order, num_samples - filter_order))
    for smpl in range(filter_order, num_samples, 1):
        # 係数保存
        coefhistory[:, smpl-filter_order] = coef
        # 予測
        pred[smpl] = np.dot(coef, data[smpl - filter_order:smpl])
        # 残差計算
        error[smpl] = data[smpl] - pred[smpl]
        # 係数更新
        coef += stepsize * error[smpl] * data[smpl - filter_order:smpl]
    return pred, coefhistory

def lms_coefhistory_experience(data, initcoef_list, stepsize):
    filter_order = initcoef_list[0].size
    listsize = len(initcoef_list)
    coefhistory = np.zeros((listsize, filter_order, data.size - filter_order))
    predhistory = np.zeros((listsize, data.size))
    for i in range(listsize):
        predhistory[i,:], coefhistory[i,:,:] = lms(data, initcoef_list[i].copy(), stepsize)
    return predhistory, coefhistory

# 係数の軌跡をプロット
def plotcoefhistory(coefhistory):
    # plt.plot(coefhistory[0,:], coefhistory[1,:],marker='o')
    plt.plot(coefhistory[0,:], coefhistory[1,:])
    
STEPSIZE = 0.05
initcoef_list = [ np.array([-1.5, -1.5]), 
                 np.array([0.5, -1.5]),
                 np.array([2.5, -1.5]), 
                 np.array([2.5,  0.5]), 
                 np.array([2.5, 2.5]),
                 np.array([0.5, 2.5]),
                 np.array([-1.5, 2.5]),
                 np.array([-1.5, 0.5]), ]

data = lineardata
predhistory, coefhistory = lms_coefhistory_experience(data, initcoef_list, STEPSIZE)
zerror, zerrorgrad = lms_experience(data, xtics, ytics)

for i in range(len(initcoef_list)):
    plotcoefhistory(coefhistory[i])
ploterror('Error', xtics, ytics, zerror)
plotgradient('Gradient', xtics, ytics, zerrorgrad)

plt.show()
In [26]:
data = sindata
_, coefhistory = lms_coefhistory_experience(data, initcoef_list, STEPSIZE)
zerror, zerrorgrad = lms_experience(data, xtics, ytics)

for i in range(len(initcoef_list)):
    plotcoefhistory(coefhistory[i])
ploterror('Error', xtics, ytics, zerror)
plotgradient('Gradient', xtics, ytics, zerrorgrad)

plt.show()
In [27]:
data = lineardata + gaussnoise
_, coefhistory = lms_coefhistory_experience(data, initcoef_list, STEPSIZE)
zerror, zerrorgrad = lms_experience(data, xtics, ytics)

for i in range(len(initcoef_list)):
    plotcoefhistory(coefhistory[i])
ploterror('Error', xtics, ytics, zerror)
plotgradient('Gradient', xtics, ytics, zerrorgrad)

plt.show()
In [28]:
data = sindata + gaussnoise
_, coefhistory = lms_coefhistory_experience(data, initcoef_list, STEPSIZE)
zerror, zerrorgrad = lms_experience(data, xtics, ytics)

for i in range(len(initcoef_list)):
    plotcoefhistory(coefhistory[i])
ploterror('Error', xtics, ytics, zerror)
plotgradient('Gradient', xtics, ytics, zerrorgrad)

plt.show()
In [29]:
data = lineardata + laplacenoise
_, coefhistory = lms_coefhistory_experience(data, initcoef_list, STEPSIZE)
zerror, zerrorgrad = lms_experience(data, xtics, ytics)

for i in range(len(initcoef_list)):
    plotcoefhistory(coefhistory[i])
ploterror('Error', xtics, ytics, zerror)
plotgradient('Gradient', xtics, ytics, zerrorgrad)

plt.show()
In [30]:
data = sindata + laplacenoise
_, coefhistory = lms_coefhistory_experience(data, initcoef_list, STEPSIZE)
zerror, zerrorgrad = lms_experience(data, xtics, ytics)

for i in range(len(initcoef_list)):
    plotcoefhistory(coefhistory[i])
ploterror('Error', xtics, ytics, zerror)
plotgradient('Gradient', xtics, ytics, zerrorgrad)

plt.show()
In [31]:
# Sined-LMS フィルター処理
def signedlms(data, initcoef, stepsize):
    num_samples = data.size
    filter_order = initcoef.size
    pred = np.zeros(num_samples)
    error = np.zeros(num_samples)
    coef = initcoef
    coefhistory = np.zeros((filter_order, num_samples - filter_order))
    for smpl in range(filter_order, num_samples, 1):
        # 係数保存
        coefhistory[:, smpl-filter_order] = coef
        # 予測
        pred[smpl] = np.dot(coef, data[smpl - filter_order:smpl])
        # 残差計算
        error[smpl] = data[smpl] - pred[smpl]
        # 係数更新
        coef += stepsize * np.sign(error[smpl]) * data[smpl - filter_order:smpl]
    return pred, coefhistory

def signedlms_coefhistory_experience(data, initcoef_list, stepsize):
    filter_order = initcoef_list[0].size
    listsize = len(initcoef_list)
    coefhistory = np.zeros((listsize, filter_order, data.size - filter_order))
    predhistory = np.zeros((listsize, data.size))
    for i in range(listsize):
        predhistory[i,:], coefhistory[i,:,:] = signedlms(data, initcoef_list[i].copy(), stepsize)
    return predhistory, coefhistory

SIGNEDLMS_STEPSIZE = 0.1
data = lineardata
_, coefhistory = signedlms_coefhistory_experience(data, initcoef_list, SIGNEDLMS_STEPSIZE)
zerror, zerrorgrad = signedlms_experience(data, xtics, ytics)

for i in range(len(initcoef_list)):
    plotcoefhistory(coefhistory[i])
ploterror('Error', xtics, ytics, zerror)
plotgradient('Gradient', xtics, ytics, zerrorgrad)

plt.show()
In [32]:
data = sindata
_, coefhistory = signedlms_coefhistory_experience(data, initcoef_list, SIGNEDLMS_STEPSIZE)
zerror, zerrorgrad = signedlms_experience(data, xtics, ytics)

for i in range(len(initcoef_list)):
    plotcoefhistory(coefhistory[i])
ploterror('Error', xtics, ytics, zerror)
plotgradient('Gradient', xtics, ytics, zerrorgrad)

plt.show()
In [33]:
data = lineardata + gaussnoise
_, coefhistory = signedlms_coefhistory_experience(data, initcoef_list, SIGNEDLMS_STEPSIZE)
zerror, zerrorgrad = signedlms_experience(data, xtics, ytics)

for i in range(len(initcoef_list)):
    plotcoefhistory(coefhistory[i])
ploterror('Error', xtics, ytics, zerror)
plotgradient('Gradient', xtics, ytics, zerrorgrad)

plt.show()
In [34]:
data = sindata + gaussnoise
_, coefhistory = signedlms_coefhistory_experience(data, initcoef_list, SIGNEDLMS_STEPSIZE)
zerror, zerrorgrad = signedlms_experience(data, xtics, ytics)

for i in range(len(initcoef_list)):
    plotcoefhistory(coefhistory[i])
ploterror('Error', xtics, ytics, zerror)
plotgradient('Gradient', xtics, ytics, zerrorgrad)

plt.show()
In [35]:
data = lineardata + laplacenoise
_, coefhistory = signedlms_coefhistory_experience(data, initcoef_list, SIGNEDLMS_STEPSIZE)
zerror, zerrorgrad = signedlms_experience(data, xtics, ytics)

for i in range(len(initcoef_list)):
    plotcoefhistory(coefhistory[i])
ploterror('Error', xtics, ytics, zerror)
plotgradient('Gradient', xtics, ytics, zerrorgrad)

plt.show()
In [36]:
data = sindata + laplacenoise
_, coefhistory = signedlms_coefhistory_experience(data, initcoef_list,  SIGNEDLMS_STEPSIZE)
zerror, zerrorgrad = signedlms_experience(data, xtics, ytics)

for i in range(len(initcoef_list)):
    plotcoefhistory(coefhistory[i])
ploterror('Error', xtics, ytics, zerror)
plotgradient('Gradient', xtics, ytics, zerrorgrad)

plt.show()

最急降下法による学習の観察

（もう色んな所でしてるけど、改めて）

  • 一旦裾野が広い方に行ってから、極値を目指しているように見える。
  • Sign-LMS は LMS と同一のステップサイズでも振動しがち。
  • Sign-LMS の方は極値付近でふらつく（LMS は極値付近では動かない）
    • Sign-LMS は勾配が消えにくいから。

→次はヘッセ行列の逆行列で勾配を修正した場合の挙動を見る

In [37]:
# Hessian-LMS フィルター処理
def hessianlms(data, initcoef, stepsize):
    num_samples = data.size
    filter_order = initcoef.size
    pred = np.zeros(num_samples)
    error = np.zeros(num_samples)
    coef = initcoef
    coefhistory = np.zeros((filter_order, num_samples - filter_order))
    S = np.eye(2)
    for smpl in range(filter_order, num_samples, 1):
        # 係数保存
        coefhistory[:, smpl-filter_order] = coef
        # 予測
        pred[smpl] = np.dot(coef, data[smpl - filter_order:smpl])
        # 残差計算
        error[smpl] = data[smpl] - pred[smpl]
        # 観測データのベクトル化
        vdata = data[smpl - filter_order:smpl].reshape((filter_order, 1))
        # 勾配
        grad = error[smpl] * vdata
        # 分散（情報）行列
        S += vdata @ vdata.T
        # ヘッセ行列の逆行列を更新
        Hinv = np.linalg.inv(S / (smpl - filter_order + 1))
        # 係数更新
        coef += stepsize * (Hinv @ grad).flatten()
    return pred, coefhistory

def hessianlms_coefhistory_experience(data, initcoef_list, stepsize):
    filter_order = initcoef_list[0].size
    listsize = len(initcoef_list)
    coefhistory = np.zeros((listsize, filter_order, data.size - filter_order))
    predhistory = np.zeros((listsize, data.size))
    for i in range(listsize):
        predhistory[i,:], coefhistory[i,:,:] = hessianlms(data, initcoef_list[i].copy(), stepsize)
    return predhistory, coefhistory
                
data = sindata
hessian_predhistory, coefhistory = hessianlms_coefhistory_experience(data, initcoef_list, STEPSIZE)
zerror, zerrorgrad, zerrorgradhessian, _ = lms_experience_hessian(data, xtics, ytics)

predhistory, _ = lms_coefhistory_experience(data, initcoef_list, STEPSIZE)

for i in range(len(initcoef_list)):
    plotcoefhistory(coefhistory[i])
ploterror('Error', xtics, ytics, zerror)
plotgradient('Gradient', xtics, ytics, zerrorgradhessian)
plt.show()

for i in range(len(initcoef_list)):
    error = predhistory[i,:] - data
    plt.plot(np.abs(error[:400]))
plt.show()

for i in range(len(initcoef_list)):
    error = hessian_predhistory[i,:] - data
    plt.plot(np.abs(error[:400]))
plt.show()

plt.hist((predhistory - data).flatten(), bins=500, histtype='step')
plt.hist((hessian_predhistory - data).flatten(), bins=500, histtype='step')
plt.show()
In [38]:
data = lineardata + gaussnoise
hessian_predhistory, coefhistory = hessianlms_coefhistory_experience(data, initcoef_list, STEPSIZE)
zerror, zerrorgrad, zerrorgradhessian, _ = lms_experience_hessian(data, xtics, ytics)

predhistory, _ = lms_coefhistory_experience(data, initcoef_list, STEPSIZE)

for i in range(len(initcoef_list)):
    plotcoefhistory(coefhistory[i])
ploterror('Error', xtics, ytics, zerror)
plotgradient('Gradient', xtics, ytics, zerrorgradhessian)
plt.show()

for i in range(len(initcoef_list)):
    error = predhistory[i,:] - data
    plt.plot(np.abs(error[:400]))
plt.show()

for i in range(len(initcoef_list)):
    error = hessian_predhistory[i,:] - data
    plt.plot(np.abs(error[:400]))
plt.show()

plt.hist((predhistory - data).flatten(), bins=500, histtype='step')
plt.hist((hessian_predhistory - data).flatten(), bins=500, histtype='step')
plt.show()
In [39]:
data = sindata + gaussnoise
hessian_predhistory, coefhistory = hessianlms_coefhistory_experience(data, initcoef_list, STEPSIZE)
zerror, zerrorgrad, zerrorgradhessian, _ = lms_experience_hessian(data, xtics, ytics)

predhistory, _ = lms_coefhistory_experience(data, initcoef_list, STEPSIZE)

for i in range(len(initcoef_list)):
    plotcoefhistory(coefhistory[i])
ploterror('Error', xtics, ytics, zerror)
plotgradient('Gradient', xtics, ytics, zerrorgradhessian)
plt.show()

for i in range(len(initcoef_list)):
    error = predhistory[i,:] - data
    plt.plot(np.abs(error[:400]))
plt.show()

for i in range(len(initcoef_list)):
    error = hessian_predhistory[i,:] - data
    plt.plot(np.abs(error[:400]))
plt.show()

plt.hist((predhistory - data).flatten(), bins=500, histtype='step')
plt.hist((hessian_predhistory - data).flatten(), bins=500, histtype='step')
plt.show()
In [40]:
data = lineardata + laplacenoise
hessian_predhistory, coefhistory = hessianlms_coefhistory_experience(data, initcoef_list, STEPSIZE)
zerror, zerrorgrad, zerrorgradhessian, _ = lms_experience_hessian(data, xtics, ytics)

predhistory, _ = lms_coefhistory_experience(data, initcoef_list, STEPSIZE)

for i in range(len(initcoef_list)):
    plotcoefhistory(coefhistory[i])
ploterror('Error', xtics, ytics, zerror)
plotgradient('Gradient', xtics, ytics, zerrorgradhessian)
plt.show()

for i in range(len(initcoef_list)):
    error = predhistory[i,:] - data
    plt.plot(np.abs(error[:400]))
plt.show()

for i in range(len(initcoef_list)):
    error = hessian_predhistory[i,:] - data
    plt.plot(np.abs(error[:400]))
plt.show()

plt.hist((predhistory - data).flatten(), bins=500, histtype='step')
plt.hist((hessian_predhistory - data).flatten(), bins=500, histtype='step')
plt.show()
In [41]:
data = sindata + laplacenoise
hessian_predhistory, coefhistory = hessianlms_coefhistory_experience(data, initcoef_list, STEPSIZE)
zerror, zerrorgrad, zerrorgradhessian, _ = lms_experience_hessian(data, xtics, ytics)

predhistory, _ = lms_coefhistory_experience(data, initcoef_list, STEPSIZE)

for i in range(len(initcoef_list)):
    plotcoefhistory(coefhistory[i])
ploterror('Error', xtics, ytics, zerror)
plotgradient('Gradient', xtics, ytics, zerrorgradhessian)
plt.show()

for i in range(len(initcoef_list)):
    error = predhistory[i,:] - data
    plt.plot(np.abs(error[:]))
plt.show()

for i in range(len(initcoef_list)):
    error = hessian_predhistory[i,:] - data
    plt.plot(np.abs(error[:]))
plt.show()

plt.hist((predhistory - data).flatten(), bins=500, histtype='step')
plt.hist((hessian_predhistory - data).flatten(), bins=500, histtype='step')
plt.show()

ヘッセ行列込みのLMSで気付いたこと

  • 解の軌跡が収束している。
  • 雑音で揺れているものの、尾根に引っかからず、極小点 (0.5, 0.5) にまっすぐ向かっているように見える。
  • こころなしか収束が早い（まだ要観察。数回の試行の平均軌跡を見れば良い？）

定式化の再確認

入力データに雑音が乗ってる。 本来の定式化では、雑音は入力と独立に乗るはず。

< 今、上で実験しているもの >

雑音 -------
　　　　　　　　　　　　　↓
入力 --> + --> < システム >  --> 観測出力
<雑音が後で入るもの> 
雑音 -------------------------
　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　    ↓
入力 --> < システム > --> + --> 観測出力

意図したものにならないと、システムは雑音に対しても係数を畳み込んでしまうから雑音と相関を持つのでは？ → 下で実験までしたけど、雑音が後に入るのは間違っていそう。観測できる入力には雑音が乗っている。観測とフィルタを畳みこむ必要がある。

In [42]:
# LMS フィルター処理（リファレンス出力に対する等化）
def lmsref(reference, data, initcoef, stepsize):
    num_samples = data.size
    filter_order = initcoef.size
    pred = np.zeros(num_samples)
    error = np.zeros(num_samples)
    coef = initcoef
    coefhistory = np.zeros((filter_order, num_samples - filter_order))
    for smpl in range(filter_order, num_samples, 1):
        # 係数保存
        coefhistory[:, smpl-filter_order] = coef
        # 予測
        pred[smpl] = np.dot(coef, data[smpl - filter_order:smpl])
        # 残差計算
        error[smpl] = reference[smpl] - pred[smpl]
        # 係数更新
        coef += stepsize * error[smpl] * data[smpl - filter_order:smpl]
    return pred, coefhistory

def lmsref_coefhistory_experience(reference, data, initcoef_list, stepsize):
    filter_order = initcoef_list[0].size
    listsize = len(initcoef_list)
    coefhistory = np.zeros((listsize, filter_order, data.size - filter_order))
    predhistory = np.zeros((listsize, data.size))
    for i in range(listsize):
        predhistory[i,:], coefhistory[i,:,:] = lmsref(reference, data, initcoef_list[i].copy(), stepsize)
    return predhistory, coefhistory

data = lineardata
reference = data
predhistory, coefhistory = lmsref_coefhistory_experience(reference, data, initcoef_list, STEPSIZE)
zerror, zerrorgrad = lms_experience(data, xtics, ytics)

for i in range(len(initcoef_list)):
    plotcoefhistory(coefhistory[i])
ploterror('Error', xtics, ytics, zerror)
plotgradient('Gradient', xtics, ytics, zerrorgrad)

plt.show()

雑音をフィルタ出力に入れた場合の考察

  • 入力に雑音が入らないので、最適解が (0.5, 0.5) ではなくなった。
  • Wikipedia 読んだら雑音をあとに入れるのはカンニングに当たる事がわかった。（つまり、実装は正しくない。ただし、等価器を事前学習する場合は話が別。）

→雑音はフィルタ入力に加えるようにして、引き続きSignedLMSの実験。

In [53]:
# Hessian-SignedLMS フィルター処理
def hessiansignedlms(data, initcoef, stepsize):
    num_samples = data.size
    filter_order = initcoef.size
    pred = np.zeros(num_samples)
    error = np.zeros(num_samples)
    coef = initcoef
    coefhistory = np.zeros((filter_order, num_samples - filter_order))
    S = np.eye(2)
    Hinv = np.eye(2)
    for smpl in range(filter_order, num_samples, 1):
        # 係数保存
        coefhistory[:, smpl-filter_order] = coef
        # 予測
        pred[smpl] = np.dot(coef, data[smpl - filter_order:smpl])
        # 残差計算
        error[smpl] = data[smpl] - pred[smpl]
        # 観測データのベクトル化
        vdata = data[smpl - filter_order:smpl].reshape((filter_order, 1))
        # 勾配
        grad = np.sign(error[smpl]) * vdata
        # 確率で重み付けした情報行列の計算
        errvar =  np.abs(error[smpl])
        # ヘッセ行列の逆行列を更新
        hgrad = Hinv @ grad 
        TAU = np.max([1 / (smpl - filter_order + 1), 0.01])
        Hinv = (1 + TAU) * Hinv - TAU * (hgrad @ hgrad.T) * laplace_distribution_underdelta(errvar, DELTA)
        coef += stepsize * hgrad.flatten()
        # S += (vdata @ vdata.T) * laplace_distribution_underdelta(errvar, DELTA)
        # Hinv = np.linalg.inv(S  / (smpl - filter_order + 1))
        # 係数更新
        # coef += stepsize * (Hinv @ grad).flatten()
    return pred, coefhistory

def hessiansignedlms_coefhistory_experience(data, initcoef_list, stepsize):
    filter_order = initcoef_list[0].size
    listsize = len(initcoef_list)
    coefhistory = np.zeros((listsize, filter_order, data.size - filter_order))
    predhistory = np.zeros((listsize, data.size))
    for i in range(listsize):
        predhistory[i,:], coefhistory[i,:,:] = hessiansignedlms(data, initcoef_list[i].copy(), stepsize)
    return predhistory, coefhistory
                
HESSIANSIGNEDLMS_STEPSIZE = SIGNEDLMS_STEPSIZE / 50
data = sindata
hessian_predhistory, coefhistory = hessiansignedlms_coefhistory_experience(data, initcoef_list, HESSIANSIGNEDLMS_STEPSIZE)
zerror, zerrorgrad, zerrorgradhessian = signedlms_experience_hessian(data, xtics, ytics)

predhistory, _ = signedlms_coefhistory_experience(data, initcoef_list, SIGNEDLMS_STEPSIZE)

for i in range(len(initcoef_list)):
    plotcoefhistory(coefhistory[i])
ploterror('Error', xtics, ytics, zerror)
plotgradient('Gradient', xtics, ytics, zerrorgradhessian)
plt.show()

for i in range(len(initcoef_list)):
    error = predhistory[i,:] - data
    plt.plot(np.abs(error[:400]))
plt.show()

for i in range(len(initcoef_list)):
    error = hessian_predhistory[i,:] - data
    plt.plot(np.abs(error[:400]))
plt.show()

plt.hist((predhistory - data).flatten(), bins=500, histtype='step')
plt.hist((hessian_predhistory - data).flatten(), bins=500, histtype='step')
plt.show()

ステップサイズを小さく取ったら発散しなくなった。

  • NLMS みたくステップサイズの収束が保証されている範囲を知りたい。
In [54]:
data = lineardata + gaussnoise
hessian_predhistory, coefhistory = hessiansignedlms_coefhistory_experience(data, initcoef_list,  HESSIANSIGNEDLMS_STEPSIZE)
zerror, zerrorgrad, zerrorgradhessian = signedlms_experience_hessian(data, xtics, ytics)

predhistory, _ = signedlms_coefhistory_experience(data, initcoef_list, SIGNEDLMS_STEPSIZE)

for i in range(len(initcoef_list)):
    plotcoefhistory(coefhistory[i])
ploterror('Error', xtics, ytics, zerror)
plotgradient('Gradient', xtics, ytics, zerrorgradhessian)
plt.show()

for i in range(len(initcoef_list)):
    error = predhistory[i,:] - data
    plt.plot(np.abs(error[:400]))
plt.show()

for i in range(len(initcoef_list)):
    error = hessian_predhistory[i,:] - data
    plt.plot(np.abs(error[:400]))
plt.show()

plt.hist((predhistory - data).flatten(), bins=500, histtype='step')
plt.hist((hessian_predhistory - data).flatten(), bins=500, histtype='step')
plt.show()
In [55]:
data = sindata + gaussnoise
hessian_predhistory, coefhistory = hessiansignedlms_coefhistory_experience(data, initcoef_list,  HESSIANSIGNEDLMS_STEPSIZE)
zerror, zerrorgrad, zerrorgradhessian = signedlms_experience_hessian(data, xtics, ytics)

predhistory, _ = signedlms_coefhistory_experience(data, initcoef_list, SIGNEDLMS_STEPSIZE)

for i in range(len(initcoef_list)):
    plotcoefhistory(coefhistory[i])
ploterror('Error', xtics, ytics, zerror)
plotgradient('Gradient', xtics, ytics, zerrorgradhessian)
plt.show()

for i in range(len(initcoef_list)):
    error = predhistory[i,:] - data
    plt.plot(np.abs(error[:400]))
plt.show()

for i in range(len(initcoef_list)):
    error = hessian_predhistory[i,:] - data
    plt.plot(np.abs(error[:400]))
plt.show()

plt.hist((predhistory - data).flatten(), bins=500, histtype='step')
plt.hist((hessian_predhistory - data).flatten(), bins=500, histtype='step')
plt.show()
In [56]:
data = lineardata + laplacenoise
hessian_predhistory, coefhistory = hessiansignedlms_coefhistory_experience(data, initcoef_list,  HESSIANSIGNEDLMS_STEPSIZE)
zerror, zerrorgrad, zerrorgradhessian = signedlms_experience_hessian(data, xtics, ytics)

predhistory, _ = signedlms_coefhistory_experience(data, initcoef_list, SIGNEDLMS_STEPSIZE)

for i in range(len(initcoef_list)):
    plotcoefhistory(coefhistory[i])
ploterror('Error', xtics, ytics, zerror)
plotgradient('Gradient', xtics, ytics, zerrorgradhessian)
plt.show()

for i in range(len(initcoef_list)):
    error = predhistory[i,:] - data
    plt.plot(np.abs(error[:400]))
plt.show()

for i in range(len(initcoef_list)):
    error = hessian_predhistory[i,:] - data
    plt.plot(np.abs(error[:400]))
plt.show()

plt.hist((predhistory - data).flatten(), bins=500, histtype='step')
plt.hist((hessian_predhistory - data).flatten(), bins=500, histtype='step')
plt.show()
In [57]:
data = sindata + laplacenoise
hessian_predhistory, coefhistory = hessiansignedlms_coefhistory_experience(data, initcoef_list,  HESSIANSIGNEDLMS_STEPSIZE)
zerror, zerrorgrad, zerrorgradhessian = signedlms_experience_hessian(data, xtics, ytics)

predhistory, _ = signedlms_coefhistory_experience(data, initcoef_list, SIGNEDLMS_STEPSIZE)

for i in range(len(initcoef_list)):
    plotcoefhistory(coefhistory[i])
ploterror('Error', xtics, ytics, zerror)
plotgradient('Gradient', xtics, ytics, zerrorgradhessian)
plt.show()

for i in range(len(initcoef_list)):
    error = predhistory[i,:] - data
    plt.plot(np.abs(error[:]))
plt.show()

for i in range(len(initcoef_list)):
    error = hessian_predhistory[i,:] - data
    plt.plot(np.abs(error[:]))
plt.show()

plt.hist((predhistory - data).flatten(), bins=500, histtype='step')
plt.hist((hessian_predhistory - data).flatten(), bins=500, histtype='step')
plt.show()

所感

  • おおむねええ感じ。さて、本当にこれが既存でないのか確かめるべきでは？
  • 有効性にあたり気になるところ
    • LMS のヘッセ行列（確率による重み付けをしない）でいいんじゃないの？→ 軽く試したら全く収束しない。要精査。
    • エントロピーや誤差は LMS より良いの？
In [ ]: